#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;

use db_parameters;
use ortho_species;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    analyse_ortholog_venn_diagram

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
my @db_param = get_db_parameters();



sub get_gene_ids_by_ortholog_type($$\@)
{
	my ($dbh, $protocol_id, $gene_id_by_type) = @_;
	my @clusters;
	for my $species($ortho_name1, $ortho_name2)
	{
		my $species_id = ($species eq $ortho_name1) ? 0 : 1;
		my $sql_cmd = <<"PL/SQLCMD";
		SELECT 
			gene_id,
			ortholog_type
		FROM 
			orthologs.ortholog_sets 
		WHERE
			ortholog_type <> 'orphan' AND
			protocol_id = $protocol_id AND
            species = '$species';
PL/SQLCMD
		for my $line (@{$dbh->selectall_arrayref($sql_cmd)})
		{
			$gene_id_by_type->[$species_id]{$line->[0]} = $line->[1];
		}
	}
}


my @gene_ids_ensp;

for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;
	next unless $protocol_name eq 'phyop';
	print STDERR "    Orthology protocol [$protocol_name]...\n";

	get_gene_ids_by_ortholog_type($dbh, $protocol_id, @gene_ids_ensp);

}



for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;
	next if $protocol_name eq 'phyop';
	print "    Orthology protocol [$protocol_name]...\n";

	# get gene_ids for this protocol
	my @gene_ids_by_type;
	get_gene_ids_by_ortholog_type($dbh, $protocol_id, @gene_ids_by_type);

	# get all gene_ids
	my @all_gene_ids;
	for my $species_id(0,1)
	{
		for my $gene_id(keys (%{$gene_ids_by_type[$species_id]}), keys (%{$gene_ids_ensp[$species_id]}))
		{
			$all_gene_ids[$species_id]{$gene_id}++;
		}
	}

	my @cnts;
	my @types = ("1 to 1", "1 to m", "m to 1", "m to m", "pseudo", "Not predicted");
	for my $species_id(0,1)
	{
		for my $type1(@types)
		{
			for my $type2(@types)
			{
				$cnts[$species_id]{$type1}{$type2} = 0;
			}
		}
	}

	# sum up counts in matrix (hash)
	print STDERR "\tSum up counts\n";
	for my $species_id(0,1)
	{
		for my $gene_id(keys %{$all_gene_ids[$species_id]})
		{
			my $type_other = $gene_ids_by_type[$species_id]{$gene_id};
			$type_other = 'Not predicted' unless (defined $type_other);
			my $type_ensp = $gene_ids_ensp[$species_id]{$gene_id};
			$type_ensp = 'Not predicted' unless (defined $type_ensp);
			$cnts[$species_id]{$type_ensp}{$type_other}++;
		}
	}
	# print out
	print STDERR "\tPrint header\n";
	for my $species_id(0,1)
	{
		my @species = ($ortho_name1, $ortho_name2);
		print "-" x 80, "\n";
		print $species[$species_id], "\n\n";
		for my $ensp_type(@types)
		{
			print "\t$ensp_type";
		}
		print "\n";

		for my $other_type(@types)
		{
			print "$other_type";
			for my $ensp_type(@types)
			{
				print "\t", $cnts[$species_id]{$ensp_type}{$other_type};
			}
			print "\n";
		}
		print "\n\n";
	}
	print "*" x 80, "\n\n";
	
}

$dbh->disconnect();

